%!TEX program = xelatex
\documentclass{article}
\setlength{\parindent}{0pt}
\usepackage{ctex}
\usepackage{tikz}
\usepackage{sectsty}
\usepackage{listings}
\usepackage{pgfplots}
\usepackage{makeidx} % 调用 makeidx 宏包，用来处理索引
\makeindex % 开启索引的收集
\usepackage{amsmath,amsthm}     
\usepackage{graphicx}     
\usepackage{hyperref} 
\usepackage{url}
\usepackage{amsfonts} 
\pgfplotsset{compat=1.16}
\usepackage{multicol}
\usepackage{cite}
\begin{document}
\section*{引言}
\subsection*{有限元方法的基本设置}
这是第一个我们实际使用有限元计算的例子。\\

我们将求解一个简单版本的泊松方程，该方程的边值为零，但右手边为非零：
\begin{align*}
-\Delta u=f & \quad\quad \mathrm in \Omega,\\
u=0 & \quad\quad \mathrm on\partial\Omega.
\end{align*}
 我们将在正方形上解这个方程 $\Omega=[-1,1]^2$您已经学习了如何在步骤1和步骤2中生成网格。在这个程序中，我们还将只考虑特定的情况$f(x)=1$，并在下一个教程程序第4步中回到如何实现更一般的情况。\\

 如果你已经学习了有限元方法的基本知识，你会记得我们需要采取的步骤，通过有限维近似来近似解$u$。具体来说，我们首先需要导出上面方程的弱形式，我们通过将方程乘以左边的测试函数$\varphi$（我们将回到下面从左边而不是从右边相乘的原因）并在域$\Omega$上积分来获得它
\begin{equation*}
    -\int_\Omega\varphi\Delta u=\int_\Omega\varphi f.
\end{equation*}
这可以分部积分:
\begin{equation*}
    \int_\Omega\nabla\varphi\cdot\nabla u-\int_{\partial\Omega}\varphi\mathbf{n}\cdot\nabla u=\int_\Omega\varphi f.
\end{equation*}
测试函数$\varphi$必须满足相同类型的边界条件（用数学术语来说：它需要来自我们寻求解的集合的切线空间），因此在边界$\varphi=0$上，我们正在寻找的弱形式读取:
\begin{equation*}
    (\nabla\varphi,\nabla u)=(\varphi,f),
\end{equation*}
其中我们使用了常用的表示法$(a,b)=\int_{\Omega}a b.$。然后，问题从适当的空间（这里是空间$H_1$）中要求函数$u$，对于该函数$u$，该语句对于所有测试函数$\varphi$都成立。\\

当然，在一般情况下，我们无法在计算机上找到这样的函数，而是寻求一个近似值$u_h(\mathbf{x})=\sum_jU_j\varphi_j(\mathbf{x})$，其中$U_j$是我们需要确定的未知展开系数（这个问题的“自由度”），$\varphi_i(x)$是我们将使用的有限元形状函数。要定义这些形状函数，我们需要以下内容：
\begin{itemize}
    \item 用于定义形状函数的网格。您已经了解了如何在步骤1和步骤2中生成和操作描述网格的对象。
    \item 一种有限元，用于描述我们要在引用单元格上使用的形状函数（在deal.II中，它始终是单位区间$[0,1]$、单位平方$[0,1]^2$或单位立方体$[0,1]^3$，具体取决于您在哪个空间维度中工作）。在第2步中，我们已经使用了$FE_Q<2>$类型的对象，它表示通过在支持点上插值来定义形状函数的常用拉格朗日元素。最简单的是$FE_Q<2>(1)$，它使用多项式次数1。在2d中，这些通常被称为双线性，因为它们在参考单元的两个坐标中的每一个坐标中都是线性的。（在1d中，它们是线性的，在3d中是三线性的；然而，在deal.II文档中，我们通常不会做出这种区分，只是简单地将这些函数称为“线性”。）
    \item 以有限元对象提供的参考单元描述为基础，枚举网格上所有自由度的DoFHandler对象。您也已经在步骤2中看到了如何做到这一点。
    \item 告诉如何从参考单元上的有限元类定义的形状函数中获取真实单元上的形状函数的映射。默认情况下，除非您明确表示不这样做，否则deal.II将为此使用（双，三）线性映射，因此在大多数情况下，您不必担心这一步。
\end{itemize}
通过这些步骤，我们现在有了一组函数$\varphi_i$，我们可以定义离散问题的弱形式：找到一个函数$u_h,i.e$，即找到上面提到的展开系数$U_j$，所以:
\begin{equation*}
    (\nabla\varphi_i,\nabla u_h)=(\varphi_i,f),\quad\quad\quad i=0\ldots N-1.
\end{equation*}
注意，我们这里遵循的约定是，所有东西都从零开始计数，这在$C$和$C++$中很常见。如果插入表示$u_h(\mathbf{x})=\sum_jU_j\varphi_j(\mathbf{x})$，并且观察到:
\begin{align*}
(\nabla\varphi_i,\nabla u_h)& =\left(\nabla\varphi_i,\nabla\left[\sum_jU_j\varphi_j\right]\right)  \\
&=\sum_j\left(\nabla\varphi_i,\nabla\left[U_j\varphi_j\right]\right) \\
&=\sum_j\left(\nabla\varphi_i,\nabla\varphi_j\right)U_j.
\end{align*}
这样，问题就变成了：找到一个向量$U$，使得:
\begin{equation*}
    AU=F
\end{equation*}
其中矩阵$A$和右手边$F$定义为:
\begin{align*}
A_{ij}&=(\nabla\varphi_i,\nabla\varphi_j), \\ 
F_i&=(\varphi_i,f).
\end{align*}
\subsection*{我们应该左乘还是右乘测试函数？}
在我们继续描述如何计算这些量之前，请注意，如果我们将原始方程从右边乘以测试函数，而不是从左边乘以，那么我们就会得到以下形式的线性系统：
\begin{equation*}
    U^TA=F^T
\end{equation*}
用一个行向量$F_T$。通过转换这个系统，这当然相当于求解:
\begin{equation*}
    A^TU=F
\end{equation*}
这在这里与上述相同，因为$A=A_T$。但通常情况并非如此，为了避免任何形式的混淆，经验表明，简单地养成从左而不是从右相乘方程的习惯（就像数学文献中经常做的那样）可以避免一类常见的错误，因为矩阵是自动校正的，在比较理论和实现时不需要转置。\\

关于本教程中的第一个例子，请参见第9步，其中我们有一个非对称双线性形式，无论我们是右乘还是左乘，它都会产生影响。
\subsection*{组装矩阵和右侧向量}
现在我们知道了我们需要什么（即：保存矩阵和向量的对象，以及计算$A_{ij}$、$F_i$的方法），我们可以看看实现这一点需要什么：
\begin{itemize}
    \item A的对象是SparseMatrix类型，而U和F的对象是Vector类型。我们将在下面的程序中看到哪些类用于求解线性系统。
    \item 我们需要一种形成积分的方法。在有限元方法中，这通常使用求积来完成，即用每个单元上一组求积点上的加权和来代替积分。也就是说，我们首先将$\Omega$上的积分分解为所有单元上的积分\begin{align*}
&A_{ij} =(\nabla\varphi_i,\nabla\varphi_j)=\sum_{K\in\mathbb{T}}\int_K\nabla\varphi_i\cdot\nabla\varphi_j,  \\
&F_{i} =(\varphi_i,f)=\sum_{K\in\mathbb{T}}\int_K\varphi_if, 
\end{align*}
然后通过正交来近似每个小区的贡献：
\begin{align*}
&A_{ij}^K =\int_K\nabla\varphi_i\cdot\nabla\varphi_j\approx\sum_q\nabla\varphi_i(\mathbf{x}_q^K)\cdot\nabla\varphi_j(\mathbf{x}_q^K)w_q^K,  \\
&F_{i}^{K} =\int_K\varphi_if\approx\sum_q\varphi_i(\mathbf{x}_q^K)f(\mathbf{x}_q^K)w_q^K, 
\end{align*}
其中$T\approx \Omega$是近似于域的三角剖分，$x^K_{q}$是单元$K$上的第$q$个正交点，$w^K_{q}$是第$q$个求积权重。做这件事需要不同的部分，我们接下来将依次讨论。
\item 首先，我们需要一种方法来描述正交点的位置$x^K_{q}$及其权重$w^K_{q}$。它们通常以与形状函数相同的方式从引用单元格进行映射，即隐式使用MappingQ1类，或者，如果您明确表示，则通过从Mapping派生的其他类之一进行映射。参考单元格上的位置和权重由从正交基类派生的对象描述。通常，人们选择一个求积公式（即一组点和权重），使得求积恰好等于矩阵中的积分；这是可以实现的，因为积分中的所有因子都是多项式，并且是通过在QGauss类中实现的高斯求积公式来实现的。
\item 然后，我们需要一些东西来帮助我们评估单元$K$上的$\varphi_i(x^K_{q})$。这就是FEValues类所做的：它使用有限元对象来描述参考单元上的$\varphi$，使用正交对象来描述正交点和权重，以及映射对象（或隐式地采用MappingQ1类），并且在位于$K$上的正交点处提供实单元$K$上的形状函数的值和导数以及积分所需的各种其他信息。
\end{itemize}
计算矩阵和右手边作为对所有单元求和(然后对正交点求和)的过程通常称为线性系统的组装，或简称组装，使用与装配线相关的单词的含义，意思是“将一组碎片，片段或元素组合在一起的行为”。\\

FEValues实际上是汇编过程中的中心类。你可以这样看待它:有限元素及其派生类描述形状函数，即无限维对象:函数在每个点上都有值。理论上我们需要这个因为我们想用函数上的积分进行分析。然而，对于计算机来说，这是一个非常困难的概念，因为它们通常只能处理有限数量的信息，因此我们通过使用在参考单元(正交对象)上定义的点映射(映射对象)到实际单元上的点来获得对正交点的求和来替换积分。从本质上讲，我们将问题简化为只需要有限数量的信息，即形状函数值和导数，正交权值，法向量等，只在有限的一组点上。FEValues类将这三个组件组合在一起，并提供关于特定单元格$K$的有限信息集。当我们组装下面的线性系统时，您将看到它的作用。\\

值得注意的是，如果您只是在应用程序中自己创建这三个对象，并自己处理信息，也可以实现所有这些目标。\\

然而，这既不会更简单(FEValues类提供了你真正需要的那种信息)，也不会更快:FEValues类被高度优化，只在每个单元格上计算你需要的特定信息;如果可以重用前一个单元格中的任何内容，那么它就会这样做，并且该类中有很多代码确保在有利的地方缓存内容。\\

本介绍的最后一部分是提到，在获得线性系统之后，使用迭代求解器对其进行求解，然后进行后处理:我们使用DataOut类创建一个输出文件，然后可以使用一种常见的可视化程序对其进行可视化。
\subsection*{求解线性系统}
对于一个有限元程序，我们最终得到的线性系统相对较小:矩阵的大小为$1089×1089$，因为我们使用的网格是$32×32$，因此网格中有$33^2=1089$个顶点。在许多后来的教程程序中，矩阵大小在数万到数十万之间的情况并不少见，并且像ASPECT这样的代码是基于deal.II构建的。我们经常用超过1亿个方程来解决问题(尽管使用并行计算机)。\\

在任何情况下，即使是这里的小系统，矩阵也比我们在本科或研究生课程中遇到的要大得多，所以问题来了，我们如何求解这样的线性系统。\\

解决线性系统的第一种典型方法是高斯消去法。这种方法的问题在于它需要的操作次数与$N^3$成正比，其中$N$是线性系统中方程或未知数的数量更具体地说，操作次数是$\frac{2}{3}N^3$，上下浮动。当$N=1089$时，这意味着我们需要进行大约8.61亿次操作。这是一个非常可行的数字，现代处理器只需不到0.1秒的时间就能完成。但很明显，这是不可能的:如果我们在线性系统中有20倍的方程(也就是说，20倍的未知数)，那么它已经需要1000--10,000秒，或者大约一个小时，将线性系统再扩大十倍。很明显，我们无法在一台计算机上再解决它。\\

我们可以通过认识到矩阵中只有相对少量的条目是非零的---也就是说，矩阵是稀疏的——来挽救这种情况。高斯消去的变化可以利用这一点，使过程大大加快;我们将在第29步中首次使用一个这样的方法---在SparseDirectUMFPACK类中实现——在接下来的几个方法中。这些高斯消去法的变化可能会使我们的问题规模达到10万或20万的量级，但超出这个量级的问题就不多了。\\

相反，我们在这里要做的是采用1952年的一个想法:共轭梯度法，或简称“CG”。CG是一种“迭代”解算器，因为它形成了收敛于精确解的向量序列;事实上，在没有舍入误差的情况下，经过$N$次这样的迭代后，如果矩阵是对称且正定的，它就会找到精确解。\\

该方法最初是作为精确求解线性系统的另一种方法而开发的，就像高斯消去法一样，但它几乎没有什么优点，几十年来基本上被遗忘了。但是，当计算机变得足够强大，可以解决高斯消去法不再有效的问题时(20世纪80年代的某个时候)，CG被重新发现，因为人们意识到它非常适合于像我们从有限元方法中得到的大型稀疏系统。这是因为
\begin{enumerate}
    \item[(i)] 它计算的向量收敛于精确的解，因此我们实际上不需要做所有的$N$次迭代来找到精确的解，只要我们对相当好的近似感到满意;
    \item[(ii)] 它只需要矩阵与向量的乘积，这对于稀疏矩阵非常有用，因为根据定义，稀疏矩阵只有$\mathcal{O}(N)$个元素，所以矩阵与向量的乘积可以用$\mathcal{O}(N)$次运算来完成，而对于密集矩阵来说，做同样的事情需要花费$N^2$个操作。因此，我们可以希望用最多$\mathcal{O}(N^2)$个操作来求解线性系统，在许多情况下，操作要少得多。
\end{enumerate}
因此，有限元代码几乎总是使用迭代求解器，如CG来求解线性系统，我们也将在本代码中这样做。(我们注意到CG方法只适用于对称的正定矩阵;对于其他方程，矩阵可能不具有这些属性，我们将不得不使用适用于更一般矩阵的迭代求解器的其他变体，如BiCGStab或GMRES。)\\

这些迭代求解器的一个重要组成部分是我们指定我们想要求解线性系统的容差---本质上，关于我们愿意在近似解中接受的误差的声明。线性系统$Ax=b$的精确解$x$的近似解$\tilde{x}$的误差被定义为$\|x-\tilde{x}\|$，但这是一个我们无法计算的量，因为我们不知道确切的解$x$。相反，我们通常认为残差，定义为$\|b-A\tilde{x}\|=\|A(x-\tilde{x})\|$，是一个可计算的测度。然后让迭代求解器计算出越来越精确的解$\tilde{x}$，直到$\|b-A\tilde{x}\|\leq\tau$。一个实际的问题是$\tau$的值应该是多少。在大多数应用中，设置:
\begin{equation*}
    \tau=10^{-6}\|b\|
\end{equation*}
是一个合理的选择。我们使$\tau$与$b$的大小(范数)成正比的事实确保了我们对解的精度的期望是相对于解的大小的。这就说得通了，如果我们让右边的$b$大10倍，那么$Ax=b$的解$x$也会大10倍，$\tilde{x}$也会;我们希望$\tilde{x}$中的精确位数与之前相同，这意味着我们也应该在残差$\|b-A\tilde{x}\|$是原始大小的十倍时终止---如果我们使$\tau$与$\|b\|$成比例，这正是我们得到的结果。\\

所有这些都将在程序中的\verb|Step3::solve()|函数中实现。正如您将看到的，用deal.II设置线性求解器非常简单。整个函数只有三行。
\subsection*{关于实现}
虽然这是使用有限元方法可以求解的最简单的方程，但该程序显示了大多数有限元程序的基本结构，并作为几乎所有后续程序都将遵循的模板。具体来说，这个程序的主类是这样的:
\begin{lstlisting}[language=C++, frame=single]
class Step3
{
  public:
    Step3 ();
    void run ();

  private:
    void make_grid ();
    void setup_system ();
    void assemble_system ();
    void solve ();
    void output_results () const;

    Triangulation<2>     triangulation;
    FE_Q<2>              fe;
    DoFHandler<2>        dof_handler;

    SparsityPattern      sparsity_pattern;
    SparseMatrix<double> system_matrix;
    Vector<double>       solution;
    Vector<double>       system_rhs;
};
\end{lstlisting}
这遵循了面向对象编程的data encapsulation，也就是说，我们尽最大努力将这个类的几乎所有内部细节隐藏在外部无法访问的私有成员中。\\

让我们从成员变量开始:它们遵循我们在上面的要点中概述的构建块，即我们需要一个Triangulation和一个DoFHandler对象，以及一个描述我们想要使用的各种形状函数的有限元对象。第二组对象与线性代数有关:系统矩阵和右侧以及解向量，以及描述矩阵稀疏模式的对象。这就是本课程所需要的(也是任何固定PDE解算器所需要的基本要素)，并且需要在整个程序中生存。与此相反，汇编所需的FEValues对象仅在整个汇编过程中才需要，因此我们在执行该操作的函数中将其创建为局部对象，并在其结束时再次销毁它。\\

其次，让我们看看成员函数。这些，也已经形成了几乎所有以下教程程序将使用的共同结构:
\begin{itemize}
    \item \verb|make_grid()|:这是一个预处理函数。顾名思义，它设置存储三角测量的对象。在后面的例子中，它还可以处理边界条件、几何图形等。
    \item \verb|setup_system()|:在这个函数中设置了解决问题所需的所有其他数据结构。特别是，它将初始化DoFHandler对象并正确地调整与线性代数有关的各种对象的大小。这个函数通常与上面的预处理函数分开，因为在一个时间相关的程序中，当网格自适应细化时，它可能至少每隔几个时间步调用一次(我们将在步骤6中看到如何做)。另一方面，在上面的预处理函数中设置网格本身只在程序开始时完成一次，因此，将其分离到自己的函数中。
    \item \verb|assemble_system()|:这是计算矩阵和右侧内容的地方，如上面的介绍中详细讨论的那样。由于处理这个线性系统在概念上和计算它的项是非常不同的，我们把它和下面的函数分开。
    \item \verb|solve()|:这就是我们计算线性系统$AU=F$的解$U$的函数。在当前的程序中，这是一个简单的任务，因为矩阵非常简单，但是只要问题不再那么简单，它就会成为程序大小的重要组成部分(例如，一旦您对库有了更多的了解，请参阅步骤20、步骤22或步骤31)。
    \item \verb|output_results()|:最后，当你计算出一个解后，你可能想用它做点什么。例如，你可能想以一种可以可视化的格式输出它，或者你可能想计算你感兴趣的量:比如，热交换器的热流，机翼的空气摩擦系数，最大桥梁荷载，或简单的数值解在某一点的值。因此，这个函数是对您的解决方案进行后处理的地方。
\end{itemize}
所有这些都由一个公共函数(而不是构造函数)，即\verb|run()|函数来完成。它是从创建该类型对象的地方调用的，并且它以适当的顺序调用所有其他函数。将此操作封装到\verb|run()|函数中，而不是从\verb|main()|调用所有其他函数，可以确保您可以更改该类中关注点分离的实现方式。例如，如果其中一个函数变得太大，您可以将其拆分为两个，并且您唯一需要关注的更改结果是在同一个类中，而不是在其他任何地方。\\

如上所述，您将在下面的许多教程程序中再次看到这种一般结构---有时函数名的拼写有变化，但本质上是功能分离的顺序。
\subsection*{关于类型的说明}
deal.II在命名空间类型中通过别名定义若干整型。(在前面的句子中，“integral”这个词被用作形容词，对应于名词“integer”。它不应该与表示曲线或曲面下的面积或体积的名词“积分”混淆。形容词“积分”广泛用于c++世界环境如“积分式”、“积分常数”,特别是等等),在这个程序中你会看到\verb|type::global_dof_index|在两个地方:一个整数类型,用于表示全球指数的自由度,也就是说,该指数中的一个特定的自由度DoFHandler对象上定义的三角(相对于特定的自由度指数在一个特定的细胞)。对于当前程序(以及几乎所有教程程序)，您将拥有数千到数百万个全局未知数(并且，对于$Q_1$元素，您将在2d和3d中分别拥有4个局部单元格和8个局部单元格)。因此，允许为全局DoF索引存储足够大的数字的数据类型是\verb|unsigned int|，因为它允许存储介于0到略大于40亿的数字(在大多数系统上，整数是32位的)。实际上，\verb|types::global_dof_index|就是这样。\\

那么，为什么不直接使用\verb|unsigned int|呢?在7.3版本之前,deal.II一直这样做。然而,deal.II 支持非常大的计算(通过步骤40中讨论的框架)，当分布在几千个处理器上时，可能有超过40亿个未知数。因此，在\verb|unsigned int|不够大的情况下，我们需要64位unsigned整型。为了实现这一点，我们引入了\verb|types::global_dof_index|，默认情况下它被定义为简单的\verb|unsigned int|，而如果有必要，可以通过在配置过程中传递一个特定的标志来将其定义为\verb|unsigned long long int|(请参阅ReadMe文件)。\\

这涵盖了技术方面。但是还有一个文档目的:在库和构建在库上的代码的任何地方，如果您看到使用数据类型\verb|types::global_dof_index|的地方，您立即知道正在引用的数量实际上是一个全局\verb|dof|索引。如果我们只使用\verb|unsigned int|(也可以是局部索引、边界指示符、材料id等)，则没有明显的意义。立即知道变量引用的内容也有助于避免错误:很明显，如果您看到\verb|types::global_dof_index|类型的对象被分配给\verb|types::subdomain_id|类型的变量，即使它们都是由无符号整数表示的，编译器也不会因此报错。\\

在更实际的术语中，这种类型的存在意味着在组装期间，我们创建一个$4×4$矩阵(在2d中，使用$Q_1$元素)我们当前所处的单元的贡献，然后我们需要将这个矩阵的元素添加到全局(系统)矩阵的适当元素中。为此，我们需要获得当前单元局部自由度的全局指标，为此我们将始终使用以下代码片段:
\begin{lstlisting}[language=C++,frame=single]
cell->get_dof_indices(local_dof_indices);
\end{lstlisting}
其中\verb|local_dof_indices|声明为
\begin{lstlisting}[language=C++,frame=single]
std::vector<types::global_dof_index> 
local_dof_indices (fe.n_dofs_per_cell());
\end{lstlisting}
这个变量的名称可能有点用词不当---它代表“当前单元格上局部定义的自由度的全局索引”---但在整个库中，保存这些信息的变量都是以这种方式通用命名的。
\subsection*{描述}
该程序解决的是二维拉普拉斯方程的边界值问题。它使用有限元方法对问题进行离散化，并利用线性有限元空间对未知函数$u$进行逼近。它将问题转化为求解一个稀疏线性系统$Ax = b$，其中$A$是系统矩阵，$b$是右侧向量，$x$是未知向量。\\

在程序的执行过程中，它会进行以下步骤：
\begin{enumerate}
    \item[1.] 创建初始网格：程序从一个初始的粗略网格开始，表示整个域Ω。

\item[2.] 网格细化：根据\verb|--refinements|参数设置的细化水平，程序将初始网格进行一定次数的全局细化。

\item[3.] 设置边界条件：根据\verb|--boundary|参数设置的类型和数值，程序为边界上的节点设置相应的边界条件。可以是迪里切特边界条件、诺依曼边界条件或混合边界条件。

\item[4.] 组装线性系统：程序根据有限元方法，将拉普拉斯方程离散化成一个稀疏线性系统$Ax = b$。它通过在每个单元上进行数值积分来计算系统矩阵$A$和右侧向量$b$的元素。

\item[5.] 求解线性系统：程序使用一个线性求解器（如共轭梯度法）来求解得到未知向量$x$，其中包含了近似的解$u$。

\item[6.] 输出结果：根据\verb|--output|参数设置的文件名，程序将数值解$u$以及其他相关信息（如网格和误差）输出到文件中，以供后续分析或可视化。
\end{enumerate}
通过调整参数，用户可以定制化问题的各个方面，例如网格细化水平、边界条件类型和数值，精确解函数等。这使得\verb|step-3|示例程序具有很好的扩展性和自定义性，适用于不同类型的二维拉普拉斯方程求解。
\subsection*{具体参数}
Number of active cells: 1024

Number of degrees of freedom: 1089

DEAL:cg::Starting value 0.121094

DEAL:cg::Convergence step 48 value $5.33692\times 10^{-13}$

\begin{figure}[h]
  \centering
  \includegraphics[width=\textwidth]{solution_3d.eps}
\end{figure}

\end{document}
